classdef FeaSelect< handle
    properties (Constant)
        
    end
    methods(Static = true)
        function kopt = GetOptimalFeature(Xcal, ycal, nwave, feaType, metal, bShow)
            switch feaType
                case 1
                    %%方便调试CARS
%                     kopt = FeaSelect.DoCars(Xcal, ycal, bShow, nwave);
                    kopt = FeaSelect.DoSPAEx(Xcal, ycal, 20, nwave, metal);
                case 2
                    kopt = FeaSelect.DoFsrfTest(Xcal, ycal, nwave, bShow);
                case 3
                    % %                     kopt = FeaSelect.DoRelieff(Xcal, ycal, bShow, ncomp);
                    % %                     kopt = FeaSelect.DoFeaRFTree(Xcal, ycal, bShow, ncomp);%%2024-07-14 改
                    kopt = FeaSelect.DoUVE(Xcal, ycal, bShow, nwave);%%2024-07-14 改
                case 4
                    kopt = FeaSelect.DoCars(Xcal, ycal, bShow, nwave);
                case 5
                    % %                     kopt = FeaSelect.DoSFS(Xcal, ycal);%%2024-07-27 取消
                    kopt = FeaSelect.DoRandomFrog(Xcal, ycal, bShow, nwave);

                case 6
                    kopt = FeaSelect.DoiRF(Xcal, ycal, bShow, nwave);
                case 7
                    kopt = FeaSelect.DoCarsNew(Xcal, ycal, bShow, nwave);                    
                case 8
                    kopt = FeaSelect.DoIriv(Xcal, ycal, bShow, nwave);%%太慢了
                otherwise
                    disp('error in DoOnePreProcessing');
                    return;
            end
            %% 2024-05-07 改
            % %             len = length(kopt{1});
            % %             kopt = kopt{1}(1: min(ncomp, len));
            
            %% 2024-05-7 改
            if length(kopt) >1
                len = length(kopt{metal});
                kopt = kopt{metal}(1: min(nwave, len));
            else
                len = length(kopt{1});
                kopt = kopt{1}(1: min(nwave, len));
            end
            
            if (len < nwave)
                disp(['len<ncomp:======len=', num2str(len), '----feaType:', num2str(feaType)])
            end
        end
        
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        function kopt = GetOptWave(XX, ycal, ins, ncomp, kkSave)
            J = size(XX, 2);
            L = length(ins);
            rmse = zeros(J, L);
            
            kopt = cell(length(ins), 1);
            for i = 1:J
                kk = kkSave{i};
                ddd = min(ncomp, length(kk));
                XXX = [ones(size(XX, 1),1), XX(:, kk(1: ddd))];
                % % % % % % % % % % % % % %                 明天这里修改
                for dd = 1 : L
                    jj = ins(dd);
                    
                    [~,~,~,~,stats] = regress(ycal(:, jj), XXX);
                    %                     eee = eee + stats(4);%% stats(4)是rmse
                    %                     eee = min(eee, stats(1));%% stats(1)是R2，因为有多个金属，rmse没法叠加，就用R2
                    rmse(i, dd) = stats(1);
                end
            end
            [M, I] = max(rmse);
            for dd = 1 : L
                kopt{dd} = kkSave{I(dd)};
            end
        end
        
        %%2024-02-08 把原DoSPA分离两个函数
        function [kkSave, XX] = GetXXSet(Xcal, K)
            x_avg = mean(Xcal);
            XX = Xcal - x_avg;
            J = size(XX, 2);
            
            kkSave = cell(J, 1);
            for i = 1 : J
                
                XX = Xcal - x_avg;
                kk = i;
                S = setdiff(1:J, kk);
                for k = 2 : K
                    xn = XX(:, kk(end));
                    
                    xn2 =  sum(xn .^2);
                    proj = xn * xn' * XX(:, S) / xn2;
                    XX(:, S) = XX(:, S) - proj;
                    Px = vecnorm(XX(:, S));
                    [M, I] = max(Px);
                    if M < 1e-20
                        break;
                    end
                    
                    kopt = S(I);
                    S = setdiff(S, kopt);
                    kk = [kk, kopt];
                end
                kkSave{i} = kk;
            end
        end
        %%SPA
        %% Xcal 24*200
        %% ycal 24*6  K < 23 ncomp = 7
        function kopt = DoSPAEx(Xcal, ycal, K, ncomp, ins)
            [kkSave, XX] = FeaSelect.GetXXSet(Xcal, K);
            kopt = FeaSelect.GetOptWave(XX, ycal, ins, ncomp, kkSave);
        end
        
        %%SPA
        %% Xcal 24*200
        %% ycal 6*24  K < 23 ncomp = 7
        function kopt = DoSPA(Xcal, ycal, K, ncomp, bShow)
            
            x_avg = mean(Xcal);
            XX = Xcal - x_avg;
            J = size(XX, 2);
            N = size(XX, 1);
            
            rmse = zeros(J, 1);
            kkSave = cell(J, 1);
            for i = 1 : J
                XX = Xcal - x_avg;
                kk = i;
                S = setdiff(1:J, kk);
                for k = 2 : K
                    xn = XX(:, kk(end));
                    
                    xn2 =  sum(xn .^2);
                    proj = xn * xn' * XX(:, S) / xn2;
                    XX(:, S) = XX(:, S) - proj;
                    Px = vecnorm(XX(:, S));
                    [M, I] = max(Px);
                    if M < 1e-20
                        break;
                    end
                    
                    kopt = S(I);
                    S = setdiff(S, kopt);
                    kk = [kk, kopt];
                end
                kkSave{i} = kk;
                
                %                 eee = 0;
                eee = 1;%%R2
                for jj = 1 : size(ycal, 2)
                    ddd = min(ncomp, length(kk));
                    XXX = [ones(size(XX, 1),1), XX(:, kk(1: ddd))];
                    
                    [~,~,~,~,stats] = regress(ycal(:, jj), XXX);
                    %                     eee = eee + stats(4);%% stats(4)是rmse
                    eee = min(eee, stats(1));%% stats(1)是R2，因为有多个金属，rmse没法叠加，就用R2
                end
                rmse(i) = eee;
            end
            
            %             [mm, I] = min(rmse);%%rmse取最小
            % %             rmse
            [mm, I] = max(rmse);%%R2去最大
            kopt = kkSave{I};
            
            if bShow == false
                return;
            end
            
            plot(1:J, Xcal);
            minn = min(Xcal, [], 'all');
            maxx = max(Xcal, [], 'all');
            hold on
            for band  = kopt
                plot([band, band], [minn, maxx], '--')
            end
            
            for i =  1: ncomp
                band  = kopt(i);
                plot([band, band], [minn, maxx], 'k-')
            end
            hold off
        end
        % % % % % % % % % % % % %
        
        
        
        
        
        %%X  32*204  Y 32*6
        function fsss = DoFsrfTest(X, y, ncomp, bShow)
            lenY = size(y, 2);
            
            fsss = cell(lenY, 1);
            for i = 1 : lenY
                yy = y(:, i);
                [idx, scores] = fsrftest(X, yy);
                
                if bShow
                    bar(scores);
                end
                fsss{i} = idx;
            end
            
            if bShow == true
                FeaSelect.plotBands(X, y, fsss, ncomp)
            end
        end
        
        % % % % % % % % % % % % % % % % % % % % % %
        %%X  32*204  Y 32*6
        function [ret, mdls] = DoFsrNCA(X, y, bShow)
            lenY = size(y, 2);
            mdls = cell(lenY, 1);
            
            for i = 1 : lenY
                yy = y(:, i);
                mdl = fsrnca(X, yy);
                mdls{i} = mdl;
                
                if bShow
                    figure()
                    plot(mdl.FeatureWeights,'ro')
                    grid on
                    xlabel('Feature index')
                    ylabel('Feature weight')
                end
            end
            
            yfit = zeros(size(y));
            for i = 1 : lenY
                yfit(:, i) = predict(mdls{i}, X);
            end
            ret = Helper.CalcAllErrors(y, yfit, bShow);
        end
        
        
        % % % % % % % % % % % % % % % % % % % % % %
        %X  32*204  Y 32*6
        function fsss = DoSFS(X, y)
            lenY = size(y, 2);
            rng('default')
            
            fsss = cell(lenY, 1);
            
            for i = 1 : lenY %%代码OK，但不用cv
                yy = y(:, i);
                c = cvpartition(yy, 'k', 10); %%6
                
                fun = @(XT,yT,Xt,yt) FeaSelect.loss(XT,yT,Xt,yt);
                %                 opts = statset('Display','iter');
                % % % % % % %                 [fs,history] = sequentialfs(FeaSelect.loss, X, yy,'cv',c,'options',opts)
                [fs,history] = sequentialfs(fun, X, yy,'cv',c)
                fsss{i} = fs;
            end
        end
        
        function ret = loss(XT, yT, Xt, yt)
            % %             mdl = fitglm(XT, yT);
            % %             ypred = predict(mdl, Xt);
            % %             ret = rms(yt - ypred);
            XX = [ones(size(XT,1),1), XT];
            Xt = [ones(size(Xt,1),1), Xt];
            beta = (XX'*XX) \ (XX'*yT);
            ypred = Xt * beta;
            ret = rms(yt - ypred);
        end
        % % % % % % % % % % % % % % % % % % % % % %
        function fsss = DoRandomFrog(X, y, bShow, ncomp)
            % % % % %             最开始的蛙跳算法，但波段比较集中，就不用了
            % % % % %             K=6; % the group number for cross validation.
            % % % % %             A=ncomp;% the maximal principle component
            % % % % %
            % % % % %             F=iRF(X, y(:, 4),10000,20,50,A, 'center')
            % % % % %             F.top10
            
            lenY = size(y, 2);
            fsss = cell(lenY, 1);
            
            % %             addpath('F:\12份土壤\matlab\libPLS_1.98')
            addpath('libPLS_1.98')
            rng('default')
            for i = 1 : lenY
                yy = y(:, i);
                
                FROG = randomfrog_pls(X, yy, ncomp);
                
                [~,indexRI]=sort(FROG.probability);
                fsss{i} = indexRI(1:ncomp);
                
                if bShow == true
                    figure
                    plot(FROG.probability)
                end
            end
            
            if bShow == true
                FeaSelect.plotBands(X, y, fsss, ncomp);
            end
            
        end
        
        % % % % % % % % % % % % % % % % % % % % % % % %
        function fsss = DoCars(X, y, bShow, ncomp)
            addpath('cars')
            
            lenY = size(y, 2);
            rng('default')
            
            
            method='center';
            A=10;
            K=6;
            N=90;%%2024-08-11 改为100
            fsss = cell(lenY, 1);
            for i = 1 : lenY
                yy = y(:, i);
                CARS = carspls2(X, yy, ncomp, K, method, N); %%原carspls与libPLS_1.98中重名，需要去除
                fsss{i} = CARS.vsel;
                
                if bShow
                    figure;
                    plotcars(CARS);
                end
            end
            
            if bShow == true
                FeaSelect.plotBands(X, y, fsss, ncomp);
            end
        end
        
        
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        %X  32*204  Y 32*6
        function fsss = DoRelieff(X, y, bShow, ncomp)
            lenY = size(y, 2);
            
            neighbors = 5;
            fsss = cell(lenY, 1);
            for i = 1 : lenY
                yy = y(:, i);
                
                [idx,weights] = relieff(X, yy, neighbors);
                if bShow == true
                    bar(weights(idx))
                    xlabel('Predictor rank')
                    ylabel('Predictor importance weight')
                end
                
                fsss{i} = idx;
            end
            
            if bShow == true
                FeaSelect.plotBands(X, y, fsss, ncomp);
            end
        end
        
        function plotBands(X, y, fsss, ncomp)
            lenY = size(y, 2);
            J = size(X, 2);
            figure
            plot(1:J, X);
            
            hold on
            minn = min(X, [], 'all');
            maxx = max(X, [], 'all');
            
            colors = [[1, 0, 0]; [0, 1, 0]; [0, 0, 1]; [0, 0, 0]; [1, 0, 1]; [0, 1, 1]];
            for i = 1 : lenY
                if ncomp < length(fsss{i})
                    bands  = fsss{i}(1 : ncomp);
                else
                    bands  = fsss{i};
                end
                for band = bands
                    plot([band, band], [minn, maxx], '-', 'color', colors(i, :))
                end
            end
            hold off
        end
        
        function plotBandsEachMetal(X, y, fsss, ncomp, setting)
            lenY = size(y, 2);
            J = size(X, 2);
            
            
            hold on
            minn = min(X, [], 'all');
            maxx = max(X, [], 'all');
            
            colors = [[1, 0, 0]; [0, 1, 0]; [0, 0, 1]; [0, 0, 0]; [1, 0, 1]; [0, 1, 1]];
            for i = 1 : lenY
                figure
                %                 plot(1:J, X);
                preType = setting(i, 1);
                XX = PreHelper.DoOnePreProcessing(X, preType, false);
                plot(1:J, XX);
                xlim([1 J])
                hold on
                if ncomp < length(fsss{i})
                    bands  = fsss{i}(1 : ncomp);
                else
                    bands  = fsss{i};
                end
                for band = bands
                    plot([band, band], [minn, maxx], '-', 'color', colors(i, :))
                end
                
                hold off
            end
            
        end
        
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        function fsss = DoLAR(X, y, bShow, ncomp)
            lenY = size(y, 2);
            fsss = {};
            
            fsss = cell(lenY, 1);
            for i = 1 : lenY
                yy = y(:, i);
                
                [b info] = lar(X, yy)
                %                 if bShow == true
                %
                %                 end
                [M. ins] = min(info.AIC);
                fsss{i} = ins;
                
            end
            
            if bShow == true
                FeaSelect.plotBands(X, y, fsss, ncomp);
            end
        end
        
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        function fsss = DoFeaRFTree(X, y, bShow, ncomp)
            lenY = size(y, 2);
            fsss = cell(lenY, 1);
            
            for i = 1 : lenY
                yy = y(:, i);
                
                tree = TreeBagger(100, X, yy, 'Method', 'regression', ...
                    'OOBPredictorImportance','On', 'MinLeafSize',3);
                
                [~, ins] = maxk(tree.OOBPermutedPredictorDeltaError, ncomp);
                fsss{i} = ins;
            end
            
            if bShow == true
                FeaSelect.plotBands(X, y, fsss, ncomp);
            end
        end
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        function fsss = DoUVE(X, y, bShow, ncomp)
            lenY = size(y, 2);
            fsss = cell(lenY, 1);
            
            %                 addpath('F:\12份土壤\matlab\libPLS_1.98')
            addpath('libPLS_1.98')
            
            for i = 1 : lenY
                yy = y(:, i);
                
                UVE = mcuvepls(X, yy, ncomp);
                RI = abs(UVE.RI);
                [RIs,indexRI]=sort(RI);
                fsss{i} = indexRI(1:ncomp);
                
                if bShow == true
                    figure
                    plot(abs(UVE.RI))
                end
            end
            
            if bShow == true
                FeaSelect.plotBands(X, y, fsss, ncomp);
            end
        end
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        function fsss = DoIriv(X, y, bShow, ncomp)
            lenY = size(y, 2);
            fsss = cell(lenY, 1);
            
            addpath('libPLS_1.98')
            
            for i = 1 : lenY
                yy = y(:, i);
                
                F = iriv(X, yy, ncomp);
                
                fsss{i} = F.SelectedVariables;
                
            end
            
            if bShow == true
                FeaSelect.plotBands(X, y, fsss, ncomp);
            end
        end
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        function fsss = DoCarsNew(X, y, bShow, ncomp)
            lenY = size(y, 2);
            fsss = cell(lenY, 1);
            rng('default')
            addpath('libPLS_1.98')
            
            for i = 1 : lenY
                yy = y(:, i);
                
                F = carspls(X, yy, ncomp, 6);
                
                fsss{i} = F.vsel;
                
            end
            
            if bShow == true
                FeaSelect.plotBands(X, y, fsss, ncomp);
            end
        end
        
         % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        function fsss = DoiRF(X, y, bShow, ncomp)
            lenY = size(y, 2);
            fsss = cell(lenY, 1);
            rng('default')
            % %             addpath('F:\12份土壤\matlab\libPLS_1.98')
            addpath('libPLS_1.98')
            for i = 1 : lenY
                yy = y(:, i);
                
                FROG = irf(X, yy, ncomp);
                
                [~, indexRI] = maxk(FROG.probability, ncomp);
                fsss{i} = indexRI;
                
                if bShow == true
                    figure
                    plot(FROG.probability)
                end
            end
            
            if bShow == true
                FeaSelect.plotBands(X, y, fsss, ncomp);
            end
        end
    end
    
end